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Abstract 

Within the framework of V-T theory of monatomic hquid dynamics, an exact equation is derived 
for a general equilibrium time correlation function. The purely vibrational contribution to such a 
function expresses the system's motion in one extended harmonic random valley. This contribution 
is analytically tractable and has no adjustable parameters. While this contribution alone dominates 
the thermodynamic properties, both vibrations and transits will make important contributions to 
time correlation functions. By way of example, the V-T formulation of time correlation functions 
is applied to the dynamic structure factor S{q,Lv). The vibrational contribution alone is shown 
to be in near perfect agreement with low-temperature molecular dynamics simulations, and a 
model simulating the transit contribution with three adjustable parameters achieves equally good 
agreement with molecular dynamics results in the liquid regime. The theory indicates that transits 
will broaden without shifting the Rayleigh and Brillouin peaks in S{q,uj), and this behavior is 
confirmed by the MD calculations. We find the vibrational contribution alone gives the location 
and much of the width of the liquid-state Brillouin peak. We also discuss this approach to liquid 
dynamics compared with potential energy landscape formalisms and mode coupling theory, drawing 
attention to the distinctive features of our approach and to some potential energy landscape results 
which support our picture of the liquid state. 

PACS numbers: 05.20.Jj, 63.50.+X, 61.20.Lc, 61.12.Bt 



I. INTRODUCTION 

Long ago, Frenkel noted that the hquid-sohd phase transition has only a small 

effect on volume, cohesive forces, and specific heat, while the liquid diffuses much more 
rapidly than the solid. From these facts he argued that the motion of a liquid atom consists 
of approximately harmonic oscillations about an equilibrium point, while the equilibrium 
point itself jumps from time to time. Following this picture, Singwi et al. 

mm 

studied a series of models in which a molecule undergoes vibrational motion for a time, 
then undergoes continuous diffusion for a time. With adjustable parameters, the model was 
able to account for incoherent parts of the neutron scattering data for water and lead. In 
discussing the theory of supercooled liquids and the glass transition, Goldstein presented 
an insightful description of the potential energy landscape, and of the vibrations about local 
minima and the transitions over barriers. From computer simulations, Stillinger and Weber 
0,0,^1 found the local minima and named them "inherent structures". They suggested 
that the equilibrium properties of liquids result from vibrational excitations within, and 
shifting equilibria between, these structures. In the same spirit, Zwanzig 1^ suggested an 
approximation for the velocity autocorrelation function, in which that function evaluated for 
the vibrations within a representative potential energy valley is multiplied by a relaxation 
factor to account for the hops between valleys. Following these ideas, an active research 
program has developed, which will be summarized in Sec. V. 

In developing V-T theory for monatomic liquids, our initial goal was to construct an 
approximate but reasonably accurate Hamiltonian, from which the partition function can 
be calculated analytically and without adjustable parameters ^|. This is the most basic 
formulation available to condensed matter theory, and it was not available for liquids prior 
to V-T theory. In constructing a tractable model for the potential energy surface, the 
following three new results were established. (1) From extensive and highly accurate analysis 
of experimental specific heat data for the elemental liquids at melt, the atomic motion was 
shown to consist almost entirely of vibrations within nearly-harmonic many-atom valleys 
[13^ . (2) From highly accurate analysis of experimental entropies of melting of the elements 



the disordering entropy was shown to be a universal constant plus small scatter 
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Q]. (3) From symmetry considerations it was concluded that the noncrystalline valleys 
belong to two classes [l^ : (a) symmetric valleys, which have remnant crystalline symmetry 



2 



and consequently have a broad range of structural potentials and vibrational frequency 
distributions, and (b) random valleys, which are of overwhelming numerical superiority 
and which all have the same structural potential and vibrational frequency distribution. 
This picture of the potential surface has since been verified by computer studies jl?! . Ifsj ] . 
Only the random valleys need to be considered further, since they completely dominate 
the potential surface. The number of random valleys is fixed by the universal disordering 
entropy of the elemental liquids |l3| . Defining an extended random valley to be the harmonic 
extension to infinity of a random valley, with intervalley intersections ignored, the leading 
order Hamiltonian is the sum of extended random valley Hamiltonians, the corresponding 
partition function is called the vibrational partition function, and any statistical mechanical 
quantity calculated using this Hamiltonian is called its vibrational contribution The 
vibrational contribution alone gives an accurate account of the thermodynamic properties of 
the elemental liquids, without adjustable parameters No other theory does this. 

Beyond the vibrational contribution, corrections for anharmonicity and for the presence of 
intervalley intersections are well defined and, though complicated, they are small [i^ . 

In V-T theory, the motion of a liquid system across the boundary between two random 
valleys is a transit. We have argued that, because of the local character of equilibrium 
fluctuations in a many particle system, each transit is accomplished by a small local group 
of atoms . Further, because of the dominance of the vibrational contribution to thermo- 
dynamic functions, the intervalley intersections must be rather sharp, and transits ne arlv 
instantaneous These predictions have since been verified by computer simulations |2l| . 
In equilibrium, transits are crucial in allowing the system to visit all the random valleys 
and thus exhibit the correct liquid entropy. Transits are the expression in V-T theory of the 
equilibrium jumps of Frenkel , and of the barrier hops of Goldstein and Stillinger 
and Weber 19 
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In practical applications, V-T theory remains accurate when anharmonicity is neglected 
entirely. Then the motion consists of the vibrational contribution, for which we have exact 
analytic equations, plus transits over perfectly sharp intervalley intersections, whose effect we 
have modeled. In this way, a one-parameter model accounts for the temperature dependence 
of the specific heat of liquids, as exemplified by the experimental data for mercury For 
a model system, equations were written for the liquid free energy, the glass free energy, 
and for the nonequilibrium rate processes which occur in the liquid-glass transition region 
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where the free energy is not defined [2^. This model quahtatively reproduces the results 
of massive rate-dependent MD calculations of Vollmayr, Kob, and Binder 2^. A one- 
parameter model for the velocity autocorrelation function gives a good account of MD 
calculations for temperatures from zero to ST^ 2M ■ Hence we have verified the applicability 
of the equations of V-T theory for the glass, for the glass-transition region, and for the liquid 
to very high temperatures. A review of V-T theory and its applications, and relations to 
other theoretical developments, has been presented 25^ . 

By now we have learned enough about V-T theory to undertake an exact formulation 
of time correlation functions. That is the purpose of the present paper. Through linear 
response theory, time correlation functions contain information on nonequilibrium processes 
(see ^], Ch. 7 and 8). Hence in contrast to thermodynamic properties, transits are expected 
to make a significant contribution to time correlation functions. Though we have previously 
studied the velocity autocorrelation function j2^, here the theory will be illustrated by the 
density autocorrelation function F{q,t), since this function probes the atomic motion on a 
more detailed level. 

In Sec. II, formally exact statistical mechanical expressions for the partition function, 
equilibrium statistical averages, and time correlation functions are derived. In Sec. Ill, we 
consider the purely vibrational contribution to equilibrium statistical averages (when transits 
are neglected) and extend the formalism to include time correlation functions. The vibra- 
tional contribution to F{q,t), and to its Fourier transform S{q,uj), are derived in Sec. IV. 
S{q, u) is of interest because it is directly observed in inelastic neutron and x-ray scatter- 
ing experiments. We then consider a model for incorporating the transit contributions to 
F{q,t) and S{q,uj), and we compare the predictions of V-T theory with MD calculations. 
Our treatment is based on classical statistical mechanics, since this is highly accurate for 
nearly all the elemental liquids, and is appropriate for comparison with MD. In Sec. V our 
conclusions, methods, and aims are related to the potential energy landscape and mode 
coupling theory programs, and we note both significant agreements and divergences. We 
summarize our results and emphasize the most important conclusions in Sec. VI. 
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II. STATISTICAL AVERAGES AND TIME CORRELATION FUNCTIONS 



For simplicity we think of an A^-atom system in a cubic box, with the motion governed by 
periodic boundary conditions. Atom K at time t has position rxit) and momentum p^it), 
for K = 1, . . . , N. The total potential energy is $({r^}), and the Hamiltonian is 

^ = E^ + ^(M)- (1) 

K 

The set {r^} represents a point in configuration space. If the potential valleys are denoted 
by the index 7, our decomposition of the potential surface means that every accessible 
configuration is in one and only one valley 7. For equilibrium statistical mechanical averages, 
only the random valleys need to be considered. 
The canonical partition function Z is 



2 = E^^' (2) 



where 7 is a random valley, and 



^■y = j dpi--- dpN j dri... drN e (3) 

The integral is over the domain of valley 7, where particle permutations are not allowed. 
The configuration integral is the same in the thermodynamic limit for each random valley, 
hence we can define the random valley partition function 2^^, as for any random valley 
7. We write the number of random valleys as , where ti^ is a parameter to be determined. 
Then 

Z = w^Zrv. (4) 

This is exact in the thermodynamic limit. The factor w'^ gives a contribution Nklnw to 
the entropy, and appears in no other thermodynamic function. Calibration to experimental 
entropy of melting for the elements gives \nw = 0.80 [13j. 

Thermodynamic properties such as energy and pressure are given by equilibrium statisti- 
cal mechanical averages of the form (A), where A{{rK, Pk}) is a time-independent quantity 
depending on atomic positions and momenta. With our decomposition of configuration 
space, the average is 

T^yldpi--- dpN dri... drN A e-^^ 
E,Idp,...dpNj^dr,...drNe-P^ ' 
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The denominator is proportional to and just as in Eqs. (2-4), every term in the in 
the denominator is equal, and the same is true of the numerator, with the result 

(A) = (6) 

where {. . .)rv is the average over the domain of a single random valley. 

We now turn to time correlation functions. Let B{t) be a quantity depending on atomic 
positions and momenta as they change in time, so that B{t) = i?({rx(t), p_ft:(t)})- For 
complex B, the corresponding autocorrelation function is 

C{t) = {B{t)B*m, (7) 

which again invoking the configuration space decomposition becomes 

E J ^Pi . . . dpN X, dn... dr^ B{t)B*iO) e"/^^ 
EJdp,...dp^J^dr,...dr^e-(^^ ' 

where the zero-time variables in B*{0) are identical with the integration variables. When 
t = 0, this quantity is a fluctuation, and is just a special case of the average considered in 
the previous paragraph, and so 

C(0) = {B{0)B*{0))r.. (9) 

When t 7^ 0, however, this is not true, because the system may be in a different random 
valley at time t than at time 0, so the variable B{t)B*{0) in Eq. (jH)) has to be evaluated along 
the equilibrium trajectory of the system. For this reason the calculation of time correlation 
functions involves difficulties that are absent when computing thermodynamic quantities. 



III. VIBRATIONAL CONTRIBUTION TO TIME CORRELATION FUNCTIONS 

In evaluating the partition function and corresponding equilibrium thermodynamic func- 
tions, we previously developed a tractable model for the vibrational contribution [isf. In 
this model, the intersections among valleys are ignored, each random valley is extended 
to infinite energy, and averages of the vibrational motion are evaluated for one extended 
random valley. The model is applied to time correlation functions in the present Section. 

For configurations in random valley 7, it is useful to write the position of atom K as 

rK{t)=RK{l) + UK{t), (10) 
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where Rk{'^) is the equihbrium position and uxit) is the displacement from equilibrium. 
The set {Rxi'j)} is the structure 7. The potential energy in valley 7 is denoted $(7), and 
is expanded about equilibrium as 

$(7) = $0(7) + $2(7) + $anh(7), (11) 

where $0(7) is the structure potential, $2(7) is the contribution quadratic in displacements, 
and ^anhil) is the remainder of $(7), the anharmonic part. The anharmonic contribution 
is almost always small, often negligible. The complicated part of $(7) is the presence on its 
boundary of intersections with neighboring valleys. But $(7) is a function in 3A^ dimensions, 
and in most directions there are no such intersections at energies accessible to the liquid. 
It is only a few directions where low lying valley- valley intersections are present, and it is 
only in these few directions where transits occur. Hence to describe the vibrational motion 
alone, we shall ignore intervalley intersections, and as a leading approximation we shall also 
neglect anharmonicity. The corresponding extended random valley has potential energy 

'^ = '^0 + ^Yl '^Ki,LjUKiULj, < UKi < OO. (12) 

KLij 

The notation is i, j = Cartesian component x, y, z, and ^xi.Lj are second position derivatives 
of $ at equilibrium. The index 7 is suppressed in Eq. (fT2| . 

The Hamiltonian for motion in an extended random valley is denoted Tivib, for vibrations, 
and is 

H„5 = '^0 + J2^ + 2Z2 ^K^,L,UK^ULr (13) 
Ki KLij 

It is useful to transform from displacements to normal modes of vibration. The normal modes 
are labeled A = 1, . . . , 3A^, their coordinates are Q\, and the transformation of displacements 
is 

UKi{t) = '^WKi,X Qxit), (14) 
A 

where for each A, WKi,\ are 3A^ real components of eigenvector A. The eigenvectors diago- 
nalize the matrix of potential energy coefficients ^Ki,Lj, and they satisfy 

^^WKi^WKi^' = Sxy, orthonormality; (15) 

Ki 

y^^WKi,xWLj,x = SxiSij, completeness. (16) 

A 
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For a given extended random valley, its microscopic geometry is encoded in its eigenvectors 
and its structure. The Hamiltonian |T3j) transforms to 



2M 2 



(17) 



where Pa = MQx, and where ux is the vibrational frequency of mode A. By definition, 
uj'l > for 3A^ — 3 modes, and = for the three modes of uniform translation. It is 
understood that the translational modes are omitted from statistical mechanical analyses. 
The structural potential $0; and the distribution g{uj) of normal mode frequencies, are each 
the same for all random valleys of a given material. 

The partition function is written in Eq. (Q. The vibrational contribution is Z^b = 
Zhrvi where Zhrv is the partition function for an extended harmonic random valley. The 
result is 

^.. = u."e-««n|^- (18) 

A 

This is subject to a correction for anharmonicity, and a correction to recover the proper do- 
main of integration for Z^, or Zrv The latter is called the boundary correction. Comparison 
with experimental data shows that both corrections are small for elemental liquids at melt. 

The time correlation function C{t) is written in Eq. (jHJ. The vibrational contribution is 
Cmb{t), obtained by replacing each random valley by its extension. The system then remains 
in a single random valley, and Eq. (jHJ reduces to 

a,5(t) = (fi(t)fi*(0)),,,, (19) 

where {■ ■ is the average for an extended harmonic random valley. The average in 
Eq. can be resolved into normal-mode time correlation functions. In quadratic order 
these functions are 

kT 

(gA(t)gA'(0)),,„ = ^'^AA' COS^At, (20) 

A 

{Pxit) P^mhrv = MkT5xX' COSOJAt, (21) 



{Qx{t)Pym,,., = 0. (22) 
An exemplary time correlation function is studied in the next Section. 
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IV. THE DENSITY AUTOCORRELATION FUNCTION 



The density autocorrelation function, or intermediate scattering function, is F{q,t), 

Fiq,t) = ^{piq,t)p{-q,0)), (23) 
where p{q, t) is the Fourier transform of the density operator, 

p(qr,t) = 5^e-^'^-^-W. (24) 

K 

First we want to evaluate the vibrational contribution to F{q,t). From Eq. (fT^ . this is 

Fv^biq, t) = ^ {piq, t)pi-q, 0)),,„ . (25) 

Because an equilibrium liquid is macroscopically isotropic, the average in (jSHl), or in 
cannot depend on the angle of q in the thermodynamic limit. But we need to treat a finite 
system, for which numerical evaluations are possible. For a finite cubic system with periodic 
boundary conditions, the q form a discrete allowed set of wavevectors consistent with the 
periodicity. We are therefore allowed to average over the directions of the allowed q. Then 
Eq. p5|l becomes 

F„t{q, t) = ^ ((p(9, t)p{-q, 0)),J^, , (26) 



where (. . . )q. is the average over the star of q. 



For motion in an extended random valley, Vxit) is written from E' 
in Eq. ()26|) works out as (the algebra is the same as for the crystal 




and the average 

m) 



KL 



KL 

where Rrl = Rk — Rl- The displacement averages simplify to 



\{[q-{uK{t)-u,mf),^^ 

WK{q) - WM + {q ■ UK{t) q ■ ^^(0)),,, . (28) 



Wxiq) is the Debye- Waller factor for atom K, 



WK{q) = \{{q-UKf)^^^. (29) 



The time dependence is entirely contained in the displacement-displacement correlation 
functions {q ■ uxit) q ■ ul{0))^^^. For a random valley, these functions vanish as t ^ oo. 
Hence it is convenient to write F{q, t) in the form 



where 



Fvib{q, t) = F^ibiq, oo) + [F^ibiq, t) - Fyib{q, oo)] , (30) 



i^.*(g,oo) = l/5^cosq-i?^ie-^-('')e-^-('')\ , (31) 

\ KL fa* 



Fvib{q,t) - F^iblq, oo) = j^l "^cos q ■ Rkl 

\ KL 



X (,-WK(Q)^-WL{q) ^^{q■^l.K{t)q■^LLiO))hrv - l] ) . (32) 

/ q* 

Notice that F^hiq, oo) is positive definite. 

The dynamic structure factor S{q, uj) is directly observed in neutron and photon scattering 
experiments, and is the Fourier transform of F{q,t), 

F(g,t)e^-*- 

oo 



S{q,u) = — F{q,t)e'-'dt. (33) 



The transform will be evaluated in the one-mode approximation, obtained by expanding 
g( -.) hrv _ 1 to leading order in Eq. (jS2I), and denoted by Fjjl^{q,t) in place of Fmb{q,t): 

F^Iiq, t) - F^ibiq, oo) = ^ / ^ cos qr ■ Rkl 

\ KL 

X e-'^-('^)e-^-('') {q ■ UK{t) q ■ t^L(O)),,, ) • (34) 



Q 



The time correlation function is transformed by Eqs. p4|l and (j^IH) to give 

{q-UKit) q-ULiO))f,,^ = 

Muj 

Then S{q, u) in the one-mode approximation is 

S^ib{q,u:) = F^ib{q,oo)6{uj) + S'^liq^uj), (36) 



Tf^q ■ '^K\ q ■ wlx cos ujxt. (35) 
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where 



A 



X e-^^('^)e-'^^('') g ■ Wkx Q-wlx) • (38) 



Finally the normal mode resolution of the Debye- Waller factor is 

2Mu 



^^(^) = L 2Mul • (^') 

A ^ 

Some physical interpretation of Syib{q,uj) is helpful. In Eq. fISBj) . FmbiQ, oo)6{uj) is the 
elastic scattering cross section at momentum transfer hq. Eq. (|H7|) expresses 5'^2('?' ^) 
sum over all modes of the one-mode cross section for scattering at momentum transfer hq, 
where fx{q)6{uj + uJx) is proportional to the cross section for scattering with annihilation of 
energy hiox in mode A, and fx{q)S{u! — ux) is proportional to the cross section for scattering 
with creation of energy hoox in mode A. The proportionality constant consists of the factors 
in front of Y2x (EH)- One observes from Eq. (jHHj) that fx{q) > 0. For most elemental 

liquids at T^T^, we expect the one-mode approximation to provide a reasonable first ap- 
proximation over the range of q where vibrational scattering is important. What is missing 
from the one-mode approximation is vibrational anharmonicity and multi-mode scattering 
events. 

rn 

We now allow for transits (the following argument is also found in [29]). When atom 
K is involved in a transit, both and change in a very short time, in such a way 
that Rx + ux remains continuous and different iable in time. A detailed model of transits 



in the atomic trajectory may be found in Chisolm et al. [2J, |2^|. Here we seek a simpler 
approximation. If the time segments between transits involving atom K are denoted jk = 
1, 2, ... , then the position of atom K at time t is RKipfKit)) + UK{'jK(t),t). F{q, t) for the 
liquid is then written, from Eq. (j^ . 

Our numerical studies provide evidence, described below in connection with Eq. (US)), that 
transits can be approximately neglected in the displacements UK{'jK{t),t). We therefore 
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make this approximation, and separately average the displacement terms in Eq. fl4Up over 
harmonic vibrations. The harmonic average can be simplified to give the result, 

Fii,{q,t) ^ 1/ J^e-^''-[«-(^-W)-«-Wl e-'^-('^)e-^-('') [1 + {q ■ q ■ 1.^(0)), + ...]) 

\ KL / q* 

(41) 

Our next step is to modify this equation so as to model the presence of transits in Rk- 

There are two ways in which transits contribute to Fiig{q,t). First, transits introduce a 
fiuctuating phase in the complex exponential in Eq. ()41|1 . and this causes additional time 
decay through decorrelation along each atomic trajectory. We model this with a relaxation 
function of the form e~"*. Second, transits give rise to inelastic scattering, in addition to 
the vibrational mode scattering already present in Eq. ()41|) . and this increases the total 
scattering cross section. We model this with a multiplicative factor. 

The leading term in Eq. PT|) gives rise to the liquid Rayleigh peak, and so is denoted 
FR{q,t). Without transits, FR{q,t) reduces to Fmb{q, oo), Eq. (p?T|l . so we model FFt{q,t) as 

Fniq, t) = C{q) F,,k{q, oo) e'^^^") *. (42) 

This function decays to zero with increasing time, in accord with the liquid property 
Fiiq{q,t) — as if: — s> oo. The relaxation rate ai{q) is expected to be around the mean 
single-atom transit rate. C{q) is positive, and greater than 1 because of the inelastic scat- 
tering associated with transits (notice the total scattering cross section is not affected by 
the factor e~"^*). 

The displacement-displacement correlation function in Eq. (PT|) gives rise to the Brillouin 
peak, and so is denoted Fsi^q, t). Without transits Fsi^q, t) reduces to Fmh{q, t) — Fyib{q, oo), 
Eq. (j34p . Empirically, we have found that this vibrational contribution alone gives an 
excellent account of the location of the Brillouin peak, and the total cross section within it, 
as compared with MD calculations and with experimental data for liquid sodium [s^. This 
suggests keeping the vibrational contribution intact, as we did in going to Eq. (jHJ, and also 
suggests that we model Fs^q^t) by 

FB{q,t) = [F„,(g,t) -F„fe(g,oo)] e-"^(«) (43) 

Note Fmhiq, t) — Fyib{q, oo) decays to zero with time, because of the decay of the vibrational 
correlation function in Eq. (jnij, and this decay gives the Brillouin peak its natural width 
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[soj l . The right side of Eq. ()43p decays faster with time, hence broadens the Brillouin peak 
from its natural width, but leaves its total cross section unchanged. 

From the above equations, our model for the dynamic structure factor is 



Siiqiq, uj) = SR{q, uj) + Ssiq, c^) 



(44) 



C{q) ai F^ibiq, oo) 



7r(a;2 + af) 



SB{q,uj) 



3kT 1 
2M 3iV 



tt2 



(uJ + UJx) + Ot2 



+ 



«2 



2 J 



(45) 



(46) 



The model has three adjustable parameters for each g, namely C(g),Q;i(g), and 02(0'). For 
a discussion of how we determine the values of these parameters, as well as other details of 
the calculation such as the addressing of possible short-time errors, see 

1.6 
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FIG. 1: S{q,uj) for motion in a single random valley at q = 0.29711 Cq^ at 24 K: solid line is 
harmonic vibrational theory (Eqs. H37II39() ). circles are MD. 



At this point let us compare theory with MD. We have done calculations for a system 
with = 500 atoms and an interatomic potential representing metallic sodium at the 
density of the liquid at melt. The potential gives an accurate account of the vibrational and 
thermodynamic properties of crystal and liquid phases, and a good account of self diffusion 
in the liquid (for summaries see 25, Im|)- When the MD system is quenched to 24 K, 
transits do not occur, and the equilibrium system moves only within a single random valley. 



13 



The circles in Fig. ^ show the corresponding S{q,u!) from MD, calculated directly from the 
defining equations (j^ . (j^ . and (jH^ at g = 0.29711 a^^ (= 0.28 Qmax where Qmax = 1.05 Cg"^ 
is the position of the first peak of the static structure factor). The theoretical equations 
(I37ll39j) for sl]l{q,uj) were evaluated for the same temperature and the same random valley, 
where the 5-functions were smoothed by making a histogram, and the results are shown 
by the solid line in Fig. ^ As the Figure shows, the theory reproduces the Brillouin peak 
from the MD results perfectly. Sl]l{q,uj) is zero at low frequencies because our system has 
no normal mode frequencies below 1.7 ps^^. From Eq. ()36|) . elastic scattering contributes 
the peak F{q, oo)6{uj) in theory and MD alike. This peak is not indicated in Fig. d To 
test the transit model, we then ran the MD system at 395 K, at which temperature the 
system has melted and is transiting frequently; the resulting S{q,u!) at the same q is shown 
by the circles in Fig. |21 We also evaluated equations ()44M6|) of the model at the same 
temperature; for details of the evaluation, including the suppression of finite- iV effects and 
the criteria used to determine the three parameters, see j22|. The individual functions 
SR{q,uj), 5*5(5,0;), and their sum Siiq{q,uj) are shown as the broken, dotted, and solid lines 
in the Figure, respectively. Just as the purely vibrational theory accounted superbly for 
the Brillouin peak at low temperatures, the inclusion of the model for transits describes the 
Rayleigh and Brillouin peaks very accurately. Notice that the location of the Brillouin peak 
has not shifted from the location predicted by the vibrational part, as predicted earlier. 
For additional calculations at different q and a detailed discussion of such matters as the 
importance of multimode scattering and the significance of the magnitudes of the model 
parameters, see 

V. COMPARISON OF THEORIES 

An extensive body of work on the dynamics of liquids has been performed in the potential 
energy landscape and mode coupling theory research programs. In this Section we discuss 
those programs and compare our aims and results with theirs. 
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FIG. 2: S{q,uj) for q = 0.29711 ^ at 395 K, from the model (solid line) and from MD (circles). 
The Rayleigh (broken line) and the Brillouin (dotted hne) contributions to the model are also 
shown separately. 

A. Potential Energy Landscape (PEL) Theories 



Instantaneous normal modes (INM) were first reported in the MD calculations of Rahman 



et al. j32]. These are the vibrational modes which diagonalize the potential energy curvature 
tensor at an arbitrary point on the potential energy surface. INM data are obtained by 
averaging over randomly chosen configurations on an equilibrium MD trajectory. Seeley and 
Keyes j33| argued that the INM averaged density of states contains information which may be 
used to construct theories of liquid state dynamics. The idea was applied to self diffusion 
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36, 
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3q by using INM data to estimate the hopping rate for Zwanzig's approximation 



12 1 . Bembenek and Laird 
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40[ classified the imaginary frequency INM as unstable modes 
(double-well modes) and shoulder modes, and correspon ding in i pro vements in the fitting of 
calculated velocity autocorrelation functions followed |4l|, |42, |42l, . Buchner, Ladanyi, and 
Stratt \^\ started with the Hamiltonian for harmonic vibrations in a system of diatomic 
molecules, and showed that the INM description is accurate only for very short times, but 
that reasonable agreement is mantained for lon ger times if the imaginary frequency modes 
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47, k8||, as well as time evolution and 



are omitted. A mean field theory was studied 

mixing of INM (4J]. Kramer et al. {5^ showed that INM data at closely spaced points 
on the MD trajectory will give accurate results for the MD time correlation functions. To 
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investigate the plane wave character of sound modes in hquid ZnCl2, Ribeiro, Wilson and 
Madden ^ carried out an INM calculation of S{q,u!). This application of INM theory 
treats instantaneous configurations as if they were equilibrium configurations, hence the 
vibrational motion used is not a solution of the equations of motion (see Keves and Fourkas 



3)- 



for water 



Additional INM calculations have been done for Rb 
58, 3, 60 



for Na 



56 



57l |. and 



In contrast to INM, quenched normal modes (QNM) are the vibrational modes evaluated 
at a local potential energy minimum (an equilibrium configuration). QNM data are usually 
obtained by averaging over the potential energy minima sampled from an equilibrium MD 
trajectory. For water at 300 K, Pohorille et al. 6^ found that the radial distribution 
calculated from QNM data are very similar to the MD results. A ppl ication of QNM theory to 
water was reviewed by Ohmine and Tanaka j63|. Cao and Voth |64.l65| argued that, because 
of the thermal fluctuations, each potential energy valley must be replaced by a self-consistent 
temperature-dependent mean field. Notice this precludes a Hamiltonian formulation of liquid 
dynamics, because a Hamiltonian requires a time and temperature independent many-body 
potential. Several authors ^ 0, discussed the issue of whether or not local information 
on the potential surface, embodied by the distribution of unstable instantaneous normal 
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modes, can be used to predict the hopping rates and barrier heights for Zwanzig's model 
of self diffusion. Rabani, Gezelter and Berne (69| introduced the cage correlation function, 
which measures the rate of change of atomic surroundings, and combined this with Zwanzig's 
approximation to account for diffusion in LJ systems. For CS2, these authors [l^l concluded 
that neither the INM nor the QNM density of states is the correct one to use in Zwanzig's 
approximation. In the language of V-T theory, the cage correlation function aims to account 
for the mean transit-induced hopping motion of a single atom. 

Currently the inherent structure (IS) research program is directed mainly toward su- 
percooled liquids and the glass transition. The standard technique is to probe the po- 
tential energy surface by quenching at random times along the equilibrium MD trajectory 
for glass forming systems. The data collected this way are equilibrium statistical averages 
and depend on the temperature of the MD calculation. Systems studied include binary LJ 
mixtures, models for water, Si02, and orthoterphenyl. The binary LJ systems apparently 
do not crystallize, but do exhibit liquid behavior above an "effective" melting tempera- 
tu. r„.. Fo. Kob.A„de.en Mn„y ...e. U fo. example. . 1.0, we,, above 
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the mode-coupling critical temperature Tmct = 0.435 [72] • The following picture emerges. 
The average inherent structure energy is virt ually constant for T > T^,, but decreases as 



T decreases below T„ 



structure energy is virtually 
6^, 3, 23, 0, 75, 76, H 



78|. Schr0der et al 



presented nu- 



merical evidence that, below a crossover temperature in the vicinity of Tmct, the system 
dynamics can be separated into vibrations around inherent structures and transitions be- 
tween them. Potential energy saddle analyses also reveal this separation of system dynamics 
T^jHjH]- Consensus is that the MD systems are in thermodynamic equilibrium at tem- 
peratures down to Tj ur-r . Structural entropy decreases toward zero as temperature decreases 



toward Tmct 



n IV 1 MCT - Ji i uci ui ai em 



. The precise character of this decrease in structural 



entropy depends on the temperature-dependent distribution of the inherent structure ener- 



gies 



B H Q, H 



The Gaussian character of this distribution at low temperatures 



has been established for binary LJ systems by numerical analysis, and has been shown on 
rather general arguments to result from the central limit theorem as — oo js^. Nu- 
merical confirmation of the Adams-Gibbs relation between the decrease in structural 
entropy and the slowing down of structural relaxation has been obtained for several systems 
fl y Q. T.e .e. .o,e of Ad^.G.bb. . .Ke IS tW. of .Ke . 



discussed by Shell and Debenedetti 



but beyond our present scope 



90 



91 



89 



Recent work on ageing of the glass is of interest 



93|. 



B. Comparison of V-T and PEL theories 

The goal of PEL theories is to describe statistical mechanical properties of liquids and 
supercooled liquids in terms of temperature-dependent statistical mechanical properties of 
the potential energy surface. INM and QNM theories work with temperature-dependent 
averages of vibrational mode data, and seek to explain time correlations such as velocity 
and density autocorrelations. IS theories work with temperature-dependent averages of 
structural potential energies and vibrational frequencies, and seek to explain such quantities 
as structural and vibrational entropy. In contrast, the goal of V-T theory is to construct 
a tractable and reasonably accurate zeroth order Hamiltonian, from which both equilib- 
rium thermodynamic properties and time correlation functions can be calculated without 
adjustable parameters, and then proceed with Hamiltonian corrections to the zeroth or- 
der theory. Ultimately this is the most basic and comprehensive formulation available to 
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theoretical physics. 

The free energy in IS theory contains the important contribution fbasin, the free energy 
of the hquid constrained to be in one of its characteristic basins. The difficulties posed by 
this constraint are examined in detail by Sastry ^76\ . and Sciortino observes that there is 
no unique way to enforce the constraint 9J|. In contrast, V-T theory incorporates this con- 
straint from the outset , in a manner consistent with all available experimental data, by 
replacing each random valley by its harmonic extension to infinity, with valley- valley intersec- 
tions ignored. This formulation naturally separates in leading order and at all temperatures 
the intravalley vibrational motion and the intervalley transit motion a separation later 
found in MD calculations at temperatures in the vicinity of Tmct (121 llSl |80|, |81|. The 
V-T formulation also separates in leading order and at all temperatures the harmonic vi- 
brational free energy from other free energy contributions, a separation again found in MD 
calculations for temperatures below melting [t^, Q • 

To test the predicted uniformity [l^ of random valleys, we calculated the structural 
potential per atom and the vibrational frequency averages (a;) , (cj^) and (ln(ci;/co'o)) for 
several monatomic random valleys with = 500 — 3000 (iTI . Il8| . The rms scatter was well 
below O.OlkTm for the structural potential per atom, and was about 0.01% of the mean 



for the vibrational frequency averages 



\. This scatter is negligible in V-T theory for 



T > Tjn. We also presumed that the scatter would vanish as ^ oo. This limiting property 
is verified for the structur al p otential per atom by the analysis of Heuer and Biichner (see 
Eq. (14) with M ^ of 

As PEL theories develop, they are increasingly concerned with anharmonicity as it af- 
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95, 



96 



97|. In the PEL 



fects the interrelations among statistical mechanical averages 
approach, anharmonicity is taken to mean all atomic motion which is not harmonic vi- 
brational motion. Anharmonicity therefore poses severe theoretical difficulties as soon as 
transits become important, i.e. at T a little above Tmct- Sciortino observes that there is 
currently no satisfactory model for anharmonicity at T above Tmct In contrast, V-T 
theory separates the correction to the zeroth order potential energy surface into two dis- 
tinct contributions ^13]: (a) the boundary correction which subtracts the potential energy 
surfaces that were extended beyond the intersections of valleys, and (b) the anharmonic 
correction which is the anharmonic potential of each random valley up to its intersection 
with a neighboring valley. This is a most useful separation because the boundary correction 
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to entropy is always negative, can be reasonably estimated, and is of increasing importance 



as temperature increases 



19|. The anharmonic free energy is quite complicated but, in 



13 



20|, just as it is in the crystal 



3. 



our experience so far, it is always small in the liquid 

Recent MD studies of complex systems exhibit behaviors characteristic of the V-T theory 
description of the potential surface. For a model of orthoterphenyl, the anharmonic entropy 
is very small at all T, the structural entropy is small and has very little temperature depen- 
dence, and the harmonic vibrational entropy is most of the total (85%) and has nearly all the 
temperature dependence (Fig. 12 of These properties are very similar to the entropy 

contribution for all monatomic elemental liquids, as expressed in V-T theory QQ. In 
addition, the excess of the mean potential energy over the harmonic vibrational contribution 
decreases for the binary LJ system at T > (Fig. 5 of 9^). This decrease is present in 
all monatomic elemental liquids at T > T^, where it has been attributed to the boundary 
contribution Q • 

That the random valleys dominate the statistical mechanics of monatomic liquids at 

n 

T > Tm has been demonstrated by MD calculations (Figs. 1 and 2 of [12|). Within the 
results of recent MD calculations for systems more complex than the monatomics, there is 
also strong evidence for the statistical mechanical dominance of macroscopically uniform 
(hence random) valleys at T > Tm- The most prominent such evidence is (a) the mean 
inherent structure energy has very little temperature dependence at each fixed density for 
T > Tm, and (b) the mean vibrational frequency distribution and its log moment are nearly 
independent of temperature at each fixed density for T > Tm- Property (a) is found in 
monatomic LJ (Fig. 2 of [t?!), in binary LJ systems (Fig. lA of ^^), Fig. 2 A lower of jzJ], 
ig. 4 top of Q, Fig. 8 of 0, Fig. 2 of 3) Fig. 5 of Fig. la of x4 Fig. la of 



93[) and in a model 
(Figs. 6 and 7 of [74 
(Figs, lb and 2c of 



or water (Fi g, l b of 



Property (b) is found in binary L J systems 
Fig. 10 of '753, Figs. 2a and 3 of j93|) and in the model for water 
) . It therefore appears that the symmmetry classification of potential 
energy valleys as random and symmetric, with their associated properties, as introduced in 
V-T theory y^J is appropriate for liquids more complex than monatomics. 
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C. Comparison of V-T and Mode Coupling theories 



Detailed summaries of mode coupling theories of liquid dynami cs ar e given by Boon 



and Yip 



99| . Hansen and McDonald [26|, and Balucani and Zoppi 



100|. Applied to the 



glass transition, mode coupling theory predicts adynamic singularity when temperature 



is lowered below the critical temperature Tj 
binary LJ system, Kob and Andersen [tiI 




101 



102j . From MD calculations for the 



104| showed that mode coupling theory is 



able to rationalize the density correlation functions at low temperatures, in the vicinity of 
Tmh t- Further tests of mode coupling theory for glassy dynamics are reviewed by Gotze 
|l05[. Mode coupling theory has been extensively applied to inelastic scattering in liquids 
at T > Tm, and this is the application for which we shall compare mode coupling and V-T 
theories. 

Mode coupling theory works with the generalized Langevin equation for F{q,t), and 
expresses the memory function in terms of processes through which density fluctuations 
decay In the viscoelastic approximation, the memory function decays with 



a g-dependent relaxation time |2a 



the combined experimental data 



persion curve in liquid Rb 
approximation, Bosse et al. 



109 
110 



99l |lOO|. This a pproxima tion provides a good fit to 



and MD data 1071 . Il08| for the Brillouin peak dis- 



(see also Fig. 9.2 of 26|). Going beyond the viscoelastic 



lllj l constructed a self-consistent theory for the longitudi- 



nal and transverse current fluctuation spectra, each expressed in terms of relaxation kernels 
approximated by decay integrals which couple the longitudinal and transverse excitations. 
This theory is in good overall agreeme nt w ith extensive neutron scattering data and MD 



calcu 



112 



ations for Ar near its triple point 



lll| . The theory was developed further by Sjogren 



113|, who separated the memory function into a binary collision part, approximated 



with a Gaussian ansatz, and a more collective tail represented by a mode coupling term. For 
liquid Rb, this theory gives an "almost quantitati ve" agreem ent with results from neutron 
scattering experiment s 1 1 f| and MD calculations [iOtI llOsj ]. More recently, inel astic x-ra.\ 
scattering measurements have been done for the light alkali metals Li [llJ ] and Na 115 . 



These data have been analyzed by mode coupling theory, and the re sulting flts to > 5'(q, a;) ar e 



115, 



117, 



118 . 



119, 



120|. 



excellent, both for the experimental data and for MD calculations 

Recent results of V-T theory exhibit the following properties for the example of liquid 
sodium at melt, (a) The vibrational contribution alone, with no adjustable parameters. 
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gives an excellent account of the Brillouin peak location at all g, and gives the natural 
width of the Brillouin peak, which is at least half the total width 3^. (b) The model 
described in Section IV, with two relaxation rates and a single multiplicative factor, can be 
made to fit S{q, uj) vs u for a wide range of q, with relaxation rates near the transit rate as 
theoretically predicted [2^. The result illustrates the important point of comparison between 
V-T and mode coupling theories: the two methods are based on different decompositions 
of the physical processes involved. While mode coupling theory analyzes F{q, t) in terms 
of processes by which density fluctuations decay, V-T theory analyzes F{q,t) in terms of 
the two contributions to the total liquid motion, vibrations and transits. These two parts 
contribute to S{q, u) as follows, (a) The extended random valley vibrational modes have 
infinite lifetimes, and represent no decay processes at all, yet they give the Brillouin peak its 
natural width, (b) In zeroth order where anharmonicity is neglected, transits are responsible 
for all genuine decay processes, and they contribute the entire width of the Rayleigh peak 
and part of the width of the Brillouin peak. 



VI. SUMMARY AND DISCUSSION 



The configuration part of an equilibrium statistical mechanical average can be thought 
of as an average over the many-particle potential energy surface. In V-T theory this average 
is constructed from integrals over individual potential energy valleys, where the entire con- 
tribution in the thermodynamic limit is from random valleys. Formally exact formulas are 
given for the canonical partition function in Eqs. © and Q, and for the time correlation 
function {B{t)B*{0)) in Eq. (jH)). For time-independent averages, which give the thermody- 
namic functions, the configuration contribution is the same for every random valley, as in 
Eq. (jni). However, the same simplification does not hold for {B{t)B*{0)), because the system 
may be in a different random valley at time t than at time 0. 

In a weighted integral over a random valley, the valley- valley intersections give only 
a small contribution, because these intersections are present in a relatively small part of 
configuration space. Hence in each such integral, only a small error is made when the actual 
random valley is replaced by its harmonic extension to infinity, and the integral is extended 
as well. We model the vibrational contribution to any equilibrium statistical mechanical 
average by that average evaluated for a system moving in one extended random valley. This 
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model has no adjustable parameters, and can be evaluated for any temperature. The neglect 
of intervalley intersections, and the incidental neglect of anharmonicity, can be corrected for 
when necessary. The vibrational contribution to the partition function is written in Eq. (|18|) , 
and to a time correlation function in Eq. (jl9p . 

By way of example, the vibrational contribution is studied for the dynamic response 
functions F{q,t) and S{q,uj). The inelastic part in the one- mode approximation, Sl]l{q,uj), 
consists of independent scattering from each of the complete set of normal vibrational modes. 
This scattering determines the location and natural width of the Brillouin peak. Vibrational 
theory is in near perfect agreement with MD calculations. Fig. 1. 

In contrast to the thermodynamic functions, time correlation functions are strongly in- 
fluenced by transits. Transits cause the system trajectory to move among random valleys, 
and in the dynamical variable B{t)B*{0), each transit contributes to the decorrelation of 
B{t) with B*{0). A simple approximation for F{q,t), Eq. ()41|). preserves this effect in the 
atomic equilibrium positions, but neglects it in the atomic displacements. This approxima- 
tion yields two results: (a) while with vibrational motion alone we have Fyib{q, oo) > 0, when 
transits are present we have Fug^q, oo) = 0, and (b) while with vibrational motion alone the 
Rayleigh peak has zero width and the Brillouin peak has its nonzero natural width, the 
effect of transits is to broaden both peaks without shifting them. These transit properties 
are incorporated into a parametrized model, Eqs. ()42M6|) . which is capable of accounting 
for MD calculations of S{q,uj) in the liquid Fig. 2. 

V-T theory rests on a new analysis of the many-body potential energy surface underlying 
the motion of a monatomic liquid, namely the symmetry classification of valleys and the 
statistical dominance and microscopic uniformity of the random valleys. The same potential 
surface is then appropriate for supercooled liquid and glass states. The glass corresponds 
to very small transit rate 2^ , vanishing on the timescale of dynamic response experiments, 
so property (a) of the preceding paragraph i mplie s that F{q, oo) is positive for the glass. 



This resu lt has been observed for real glasses 
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12ll | and also for computer models of glasses 



1221 . From V-T theory, the value of this positive constant is Fyih{q, oo) for an extended 



random valley. It has also been shown p reviously that a vibrational analysis agrees with 



MD calculations for S{q,uj) in a LJ glass 12J, |l25[. The similar result for sodium is shown 



here in Fig. 1, where the vibrational theory is specifically for an extended random valley. 
Then, according to V-T theory, this vibrational contribution has application to the liquid 
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state, and leads to the result shown in Fig. 2. 

Finally, we have also compared and contrasted our results with those of the potential 
energy landscape (PEL) and mode coupling theory (MCT) approaches to liquid dynamics. 
The various types of PEL theory seek to describe the statistical mechanics of liquids in terms 
of temperature-dependent properties of their potential energy surfaces, while MCT describes 
the liquid in terms of different treatments of the Langevin equation for F{q,t). In contrast, 
our goal is to produce a theory with an accurate and tractable zeroth order Hamiltonian, 
from which equilibrium and nonequilibrium quantities can be computed without adjustable 
parameters, and then make corrections to this Hamiltonian. In so doing, we have predicted 
and verified several significant properties of the potential energy surface, such as the classifi- 
cation of valleys described in the previous paragraph, making the analysis of liquid dynamics 
considerably simpler. We have then identified two distinct corrections to the zeroth order 
Hamiltonian, namely the boundary correction which accounts for valley-valley intersections, 
and the anharmonicity of the random valleys. Of these two, the boundary correction is by 
far the most important, and is not too difficult to estimate. This helps address a significant 
problem faced recently by PEL theories. Further, some studies in the PEL program are 
revealing properties of the potential energy surface predicted by V-T theory independently 
of our efforts. In contrasting our treatment of F{q, t) with MCT, we noted that our theo- 
ries are based on different decompositions of the relevant physical processes: MCT analyzes 
F{q,t) in terms of the decay of density fluctuations, while V-T analysis is in terms of the 
two contributions to the actual motion of particles in a liquid, vibrations and transits. We 
argue that the direct connection of this analysis to the motion of the particles, which is of 
universal applicability to all equilibrium and nonequilibrium quantities, provides a basis for 
a unified treatment of all processes in the liquid and supercooled states. 

There is much still to learn about the physical nature of transits. We have observed from 
the outset that transits must be "correlation controlled," and not just thermally activated 
jl^l . This is because, as mentioned above, the valley- valley intersections are present in a 
relatively small part of configuration space, so that the transit rate is limited by the time 
required for the system to find a transit pathway. The highly correlated motion of the 
atoms involved in a transit is revealed in MD calculations [21]. Hence transits provide a 
new and challenging problem in the motion of many particle systems, and the solution of 
this problem will have a key role in the nonequilibrium properties of liquids, through their 
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time correlation functions. 
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